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Abstract 

This paper addresses the problem of the determination of the subcritical aeroelastic 
response and flutter instability of nonlinear two-dimensional lifting surfaces in an 
incompressible flow-field via indicial functions and Volterra series approach. The related 
aeroelastic governing equations are based upon the inclusion of structural and damping 
nonlinearities in plunging and pitching, of the linear unsteady aerodynamics and 
consideration of an arbitrary time-dependent external pressure pulse. Unsteady aeroelastic 
nonlinear kernels are determined, and based on these, frequency and time histories of the 
subcritical aeroelastic response are obtained, and in this context the influence of the 
considered nonlinearities is emphasized. Conclusions and results displaying the implications 
of the considered effects are supplied. 


Nomenclature 


a Dimensionless elastic axis position measured from the mid-chord, positive aft 
c Chord length of 2-D lifting surface, 2b 

ccK,K Damping and stiffness coefficients in plunging and pitching (i=l,2,3 - linear, 

quadratic, cubic), respectively 
C La Lift-curve slope 

C(k\ F{k\G{k) Theodorsen’s function and its real and imaginary counterparts, respectively 
h Plunging displacement and its dimensionless counterpart, {h/b), respectively 
h H n-th order Volterra kernel in time, and its Laplace transformed counterpart, respectively 

n 1 n 

I a ’ r a Mass moment of inertia per unit wingspan and the dimensionless radius of gyration, 

(/„ /mb 2 )' 2 , respectively 

l a ,m a Dimensionless aerodynamic lift and moment, (L a blmUl) and (m a b 2 / / a Ul), 
respectively 

L , M Total lift and moment per unit span 

a 1 a 

L b ,l b Overpressure signature of the N-wave shock pulse and its dimensionless counterpart, 

( L h b/mUl ), respectively 

m , H Airfoil mass per unit length and reduced mass ratio, (m/npb 2 ), respectively 
N Load factor, l + h'/g 

Pm ’Pm Peak reflected pressure amplitude and its dimensionless counterpart, (P m b/ milt), 
respectively 

r Shock pulse length factor 

Sj , X Laplace transform variable and Laplace operator, respectively, s } = ikj ; t ‘ = - 1 
S y Static unbalance about the elastic axis and its dimensionless counterpart, S a /mb, 
respectively 

t ,r 0 ,rTime variables and dimensionless counterpart, (Uj/b), respectively 

t p ,T p Positive phase duration, measured from the time of the arrival of the pulse, and its 

dimensionless value, respectively 
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TF Transfer function 

, y Freestream speed and its dimensionless counterpart, (U„/b(t) a ) 
x(t) Time-dependent external pulse (traveling gusts and wake blast waves) 
y(r) Response of the considered degree of freedom (pitch a and/or plunge h) 
a Twist angle about the pitch axis 

£ h ,£ a Structural damping ratios in plunging (c h /2m0 ) h ), and pitching (< c a /2I a 0 ) a ), 

respectively 
p Air density 

<p(r) Wagner’s indicial function 

CD , k Circular and reduced frequencies, (cob IUj), respectively 
(O h ,(O a Uncoupled frequencies in plunging and pitching, ( K h /m ) and ( K a jl a ) , 

respectively 

To Plunging-pitching frequency ratio, {(0 h / (O a ) 

Superscript 

( " ) Variables in Laplace transformed space 

(•),(• ) Derivatives with respect to the time t , and the dimensionless timer , respectively 



I. Introduction 

It is a well-known fact that within the linearized approach of the aeroelasticity discipline, it is 
possible to obtain the divergence and the flutter instability boundaries, and also to get the 
linearized subcritical aeroelastic response to time-dependent external pulses, such as blast and 
sonic boom pressure signatures and gust loads. On the other hand, the nonlinear approach of the 
aeroelastic problem can provide important information: i) on the influence of the considered 
nonlinearities on the subcritical aeroelastic response, and ii) about the character of the instability 
flutter boundary, i.e. benign or catastrophic one. In other words, such an approach gives the 
possibility of determining in what conditions the flutter speed can be exceeded without the 
occurrence of a catastrophic failure of the wing, in which case the flutter is benign, as well as the 
conditions in which undamped oscillations may appear at velocities below the flutter velocity, in 
which case the flutter is catastrophic. In addition, the considered nonlinearities play a great role 
on the subcritical aeroelastic response of lifting surfaces. Due to the strong implications of 
various nonlinearities on the highly flexible lifting surfaces, their related aeroelastic phenomena 
cannot longer be analyzed only within the standard linearized aeroelasticity theory. Aircraft wing 
structures often contain nonlinearities, which affect their aeroelastic behavior and performance 
characteristics, and flutter boundaries. For these reasons, in order to investigate the aeroelastic 
behavior of lifting surfaces in the vicinity of the flutter boundary, the aeroelastic governing 
equations have to be necessarily considered in nonlinear form. 

This investigation concerns the aeroelastic response of two-dimensional nonlinear lifting 
surfaces exposed to an incompressible flow field and subjected to an external pressure pulse 
Based on Volterra’s functional series approach 4 ' 8 pertinent information about the effects of 
nonlinearities on either the aeroelastic response in the subcritical flight speed regime, and their 
implication on the flutter boundary are supplied. 

The advantage of the technique based on Volterra’s series and indicial function (Lomax 9 , 
Bisplinghoff 10 , Marzocca et al. n ) consists, among others, on the possibility to investigate, within 
a rigorous theoretical basis, the aeroelastic systems featuring a wide class of nonlinearities. As a 
limiting case, based upon the first order Volterra kernel the study of the linear aeroelastic 
stability of the systems can be carried out. This methodology can encompass the case of an 
arbitrary number of degrees of freedom and at the same time is conceptually clearer, 
computationally simpler and can provide more accurate and realistic results as compared to the 
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conventional techniques used in nonlinear aeroelastic systems based on perturbation and multiple 
scale methods. In addition, this method, as shown in this paper, features a faster convergence as 
compared to the other methodologies. Moreover, this method does not experience the limitations, 
usually characterizing the other techniques, such as the Hilbert transform 12 ' 13 , developed to 
identify nonlinear systems from the first order frequency response function (FRF), or the phase 
plane methods that can describe the motion as just a function of two variables. In contrast to 
these methods that are suitable mainly to one degrees of freedom systems, Volterra series 
approach overcomes the shortcomings facing the other methods, consisting of the impossibility 
to cope with exchange of energy between the different mode frequencies. 

Toward the end of determining the nonlinear unsteady aeroelastic kernels, the harmonic probing 
algorithm, referred to as the method of growing exponentials advanced by Bedrosian and 
Rice 14,15 , and the multidimensional Laplace transform 16 will be used. 

In addition to the aeroelastic response and determination of the flutter instability boundary, 
Volterra Series considered in conjunction with this nonlinear aeroelastic model can be used to 
study the conditions rendering the flutter boundary a benign or a catastrophic one (Librescu ' ). 
Moreover, when the closed-loop dynamic response of actively controlled lifting surface is 
analyzed, also the feedback control forces and moments should be included (Librescu 1718 , Gem 
and Librescu 19 , Librescu and Na 20 , Van Trees 21 , Chua and Ng 22 ) and Volterra’s series approach 
can still be applied towards the control purpose. 

Volterra’s series approach provides a firm basis for the treatment of the nonlinear subcritical 
aeroelastic response, in the sense that it supplies an explicit relationship between the input (any 
type of time-dependent external pulses, i.e. blast load, sonic-boom, gust loads) and its response. 
With the so-called Volterra Kernel identification scheme, the modeling of an aeroelastic system 
using this approach becomes feasible. However, this methodology requires determination for 
each specific flight conditions of the corresponding nonlinear kernel of the Volterra’s series. For 
this reason, in order to define the appropriate aerodynamic loads, the recent interest in the 
modeling of unsteady nonlinear aerodynamics by this approach has been focused on the 
identification of Volterra’s kernels in the time domain (Silva 23-26 ), and in the frequency domain 
(Tromp and Jenkins 27 ). 
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A number of fundamental contributions related to Volterra’s series, developed by outstanding 
mathematicians (Volterra 4 , Wiener 5 ) and used mainly in electrical engineering 6 ' 8 , are already 
available. 

The original studies on functional series by Volterra , have been continued in the works by 
Volterra himself, and of those by Rugh 6 , Schetzen 7 , Boyd 8 . These concepts have been mainly 
used in the general nonlinear system theory. 

Originally, the method of Volterra series and Volterra kernel identification were developed to 
identify the nonlinear behavior in electrical circuits. In the aerospace field, the fundamental 
contributions were brought by Silva -5 , who has shown that the method is also applicable to 
aeroelastic systems (aerodynamic reactions and forced structural model). Silva’s pioneering 
work 23 ' 26 in this area has opened a very promising way of modeling and approaching nonlinear 
aeroelastic systems. 

II. Basic Concepts 

Having in view that for nonlinear systems the superposition principle is not applicable, and 
having in view the different types of responses induced by unsteady aerodynamic loads and 
external excitation, a combination of transfer functions is used. For the nonlinear aeroelastic 
systems these transfer functions and the time-histories response in time and frequency domains 
are determined by taking the multi-dimensional Laplace transform of Volterra kernels of the 
related aeroelastic system via a Mathematica® code developed by these authors -8 . Our approach 
intended to address the subcritical response of the nonlinear aeroelastic governing equations, is 
based on its exact representation as an infinite sum of multidimensional convolution integrals, 
the first one, (i.e. the linear kernel) being the analogous to the linear indicial aeroelastic function. 
The full nonlinear aeroelastic response will be composed of additional higher-order 
contributions. In the frequency domain, if the nonlinear function governing a system is 
“smooth”, then for small inputs the system must be asymptotically linear 6 . One of the key issues 
is to determine, corresponding to the considered type of structural, damping and aerodynamic 
nonlinearities, the pertinent Volterra's kernels. When also the active control is implemented, the 
corresponding Volterra’s kernel should also be determined. 
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III. The Theory 

In an attempt to make the paper reasonably self-contained, several elements associated with the 
indicial functions and Volterra’s series as applied to aeroelastic system, will be supplied here. 

A. Indicial Theory and Aerodynamic Loads 

Using the aerodynamic indicial functions corresponding to the transient aerodynamic reaction 
to a step pulse, the aerodynamic forces and moments induced in any maneuver and flight regime 
can be determined. Aerodynamic forces and moments acting on a rapidly maneuvering aircraft 
are, in general, nonlinear functions of the motion variables, their time rate of change, and the 
history of the maneuvering (Tobak and Chapman' 9 ). However, in this study, the linear 
aerodynamic theory is adopted. Once the response of the system to a step change in one of the 
disturbing variables (i.e. the indicial response) is known, the indicial method permits 
determination of the response to an arbitrary schedule of disturbances. There is a critical value of 
the flight speed, referred to as the flutter speed, above which the steady motion becomes 
unstable. In a nonlinear aeroelastic system the flutter phenomenon corresponds to the instability 
known as the Hopf bifurcation, resulting in the case of supercritical Hopf bifurcation in finite 
amplitude oscillations, and in the case of the subcritical Hopf bifurcation in oscillations with 
increasing amplitudes, even if the system operates before reaching the flutter speed 30 3 '. 

We need to mention that within the nonlinear indicial theory 33 , the response of a nonlinear 
system to an arbitrary input can be constructed by integrating a nonlinear functional that involves 
the knowledge of the time-dependent input and the kernel response. Whereas, within the linear 
indicial theory the linear kernel or linear impulse response can be convoluted with the input to 
predict the output of a linear system, the nonlinear indicial theory constitutes a generalization of 
this concept. It can also be shown that the traditional Volterra- Wiener theory of nonlinear 
systems constitutes a subset of the nonlinear indicial theory. It should also be mentioned that the 
nonlinear unsteady aerodynamics valid throughout the subsonic incompressible/compressible, 
transonic and supersonic flight speed regimes can be used and determined via the use of 
nonlinear indicial functions in conjunction with the Volterra’s series approach. 

B. Volterra Functional Series Theory 

As it was shown (Rugh 6 , Schetzen 7 ), within Volterra's series approach the full response in the 
time domain, y(t), of the nonlinear systems with memory can be cast as: 


7 


( 1 ) 


y(0=Ly*(0* 


Jt= o 


where, y k (/) is expressed as: 

y*(0= JjJ£^fr-T If i-r a ,..i-Tjn J <T i yr l . 

Jt /imu 


By a change of variables, it is possible to express Eq. (2) in contracted form as: 


y*(0= 

/(: fiinw 


)Y[x{t-x i YT r 


(2) 


(3) 


It is assumed that *(/)= 0 for r < 0 , implying that the system is causal. 

With this restriction, all the integrals in the subsequent discussions are different from zero over 
the time range [0,°o). Restricting the development of Eq. (3) to the first three terms one obtains: 


y{t ) = J h x (t, )x(r - r, )dx x + JJ h 2 (r, , r 2 Mr - r, M t - t 2 )dT x dr ? 

On the other hand, the response of the system can be expressed also in the frequency domain. 
Volterra’ series is essentially a polynomial approximation of the system, extension of Taylor 
series to systems with memory, while Volterra’ s kernels /i,(M a direct extension of the 
impulse response concept of the linear system theory to multiple dimensions (Volterra 5 , Rugh 6 , 
Schetzen 7 , Boyd 8 ). Consequently, a multidimensional analogue of the impulse response can be 
used to characterize a nonlinear system (Silva -3 -6 ). 

Having in view that the memory of aeroelastic systems is not infinite and, at the same time, the 
time-dependent external excitations, such as impulse, gust, blast and sonic-boom pressure 
signatures are non persistent, (in the sense that their effect decay as time unfolds), it is possible 
to characterize a nonlinear aeroelastic system via Volterra series. This fact is reflected in the 
interpretation of the Volterra kernels as higher order impulse response functions, i.e. 


h( T,,---,T n )->0 as -4' 


We will use the definition of the nonlinear transfer function or higher-order impulse response 
functions namely: 




( 5 ) 


as well as of their inverted counterparts: 
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Once Volterra’s kernels are known, the response of the nonlinear aeroelastic system can fully be 
identified. As demonstrated in the Schetzen works 7 , without loss of generality, the kernels will 
be taken as symmetric, in the sense of H n (s l ,s 2 ,...)= and a similar relationship 

being valid for h n , as well. 


If we focalize the attention on the linear system, the Laplace transform £ of the first term of Eq. 
(4) yields the familiar Laplace domain expression T(s) = //(.s)x(.s), where Y (s\ H (s), X (s) are 
the Laplace transforms of y(x), h{r), x{r) , respectively. H{s ) is the transfer function TF of the 
system. It is well known fact that for the linear system, either the first transfer function or the 
first kernel in time /i(t) encode all the information about the aeroelastic system. Moreover, as is 
well known, if the system is linear, i.e. superposition principle holds valid and is time invariant, 
the external load is uniquely related to the response by a convolution integral. With the use of 
functional series, i.e. the Volterra series, this functional representation can be extended to 
nonlinear systems. The comparison between the prediction of the linear aeroelastic responses of 
2-D lifting surface in incompressible flow field based on the Volterra’s series approach (using 
Theodorsen’s function C{k )) and on the exact solution, based on convolution integrals (using 

Wagner’s function </>(r)) is presented in Fig. 1. 

The excellent agreement of these two predictions, assess both the accuracy of the aeroelastic 
model and also the power of the methodology that combines Volterra’s series and the indicial 
function. 

IV. Mathematical Formulation 

A. General Theory for Wing Sections Including Structural and Damping 

Nonlinearities in Plunging and Pitching 

The aeroelastic governing equation of motion for 1- and 2-DOF lifting surfaces featuring 
structural and damping nonlinearities, that include the stiffnesses and the damping in plunging 
and pitching, will be analyzed next. In this sense: a 1-DOF lifting surface featuring purely 
plunging motion and a 2-DOF lifting surface featuring inertial and aerodynamic coupling in 
plunging h and pitching cc will be analyzed. As previously mentioned, the unsteady aerodynamic 
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is considered linear. A harmonic time dependent external concentrated load is also included in 
the analysis. This can be considered to correspond, for example, to an engine mounted on an 
aircraft wing. As a resuit, a harmonic type loading due to the engine oscillations has to involve 

the motion of the wing. 

As a characteristic of this approach, the transfer functions of the system would exist and would 
be the same for any excitations 34 ' 36 , such as impulse, gusts, airblast or sonic-booms (random or 
deterministic ones). This is due to the fact that transfer functions are independent of the input to 
the system, being a characteristic of the system itself. As a reminder, the validity of this method 
is based on the use of continuous polynomial type nonlinearities. For nonlinear ordinary 
differential systems, there are in general, an infinite number of Volterra kernels. In practice, one 
can handle only a finite number of terms in the series, which leads to the problem of truncation 
accuracy. However, Wiener 5 suggest that the first terms of the series may be sufficient to 
represent the output of a nonlinear system if the nonlinearities are not too strong. 

The use of the multidimensional Laplace transform as a function of several variables is a tool 
useful in stationary nonlinear system theory. The multivariable convolutions can be represented 
in terms of products of Laplace transforms. 

It is well known that the nonlinear aeroelastic systems cannot be described by a simple transfer 
function for two main reasons: a) the response encode both the unsteady aerodynamic loads and 
the external excitation effects, and b), in the nonlinear case the superposition principle is not 
applicable. It is also well known that any time-dependent external excitation, i.e. periodic or 
otherwise, can be represented, to an arbitrary degree of accuracy, by a sum of sinusoidal waves . 
In this context, if the external load is expressed in terms of multiple sinusoidal forms (for 
example, a traveling gust loads) this is easily convertible in the exponential form, i.e.: 

u{t)=Acos{a> A t)+Bcos(o} B t)^u{t)=^(e SA '+e^‘)+^(e SB '+e v ), ( 7 ) 

where s A = ico A and s B = ico B . For clarity of exposition, it is convenient to adopt this approach 
for a system with one degree of freedom (1-DOF). These results have more general bearing and 
can be extended for systems with multi-degrees of freedom (M-DOF). In fact, by using the 
classical approach of the one dimensional frequency response function, it is possible to derive an 
analytical form of the multi-dimensional frequency response characteristics of nonlinear systems. 
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The systems based on 1-DOF (plunging h) and 2-DOF (pitching a and plunging h), will be 
considered in the next sections. 

B. Plunging Airfoil Motion in an Incompressible Flow Field 

The nonlinear equation of an airfoil featuring plunging motion can be expressed as: 

+u*wr )-£.m= *,(»). (« 

i=l 

where i defines the degree of the considered nonlinearity. In the numerical simulations, i will 
assume the values 1,2,3, implying linear, quadratic and cubic nonlinearities. In addition, m is 
mass parameter and c; n , k> u , are the damping and stiffness parameters, respectively, that are 

associated to the damping and deflection in plunging corresponding to the i-th power. In the right 

hand side member of these equations, L b { t) denotes the external time-dependent load acting on 
the rigid wing counterpart. In Eq. (8) the related unsteady aerodynamic lift is represented as a 
function of the plunging degree of freedom h, only: 

L. (t ).-C u pul £_*( t - t 0 yjT' -ipc^ulk’.. (9) 

The non-circulatory component present in Eq. (9) has been represented in term of convolution 
integral of the indicial Wagner’s function. 

In order to explain how this methodology works, let us to determine, in terms of Volterra series, 
how a system responds to a harmonic or periodic time-dependent load. Let consider a periodic 
external excitation in the form: 

a)=iv"- <10) 

j=i 

As is well known, the information acquired by the case of the response to a harmonically time- 
dependent load can be used to obtain the response to any time-dependent external excitation. In 
fact, considering the case of a concentrated load arbitrarily located in the x, y plane of the wing, 
we have: 

u{x,y,t)=AS{x-x 0 ,y-y 0 )e‘ o *, (11) 

where <5(), x 0 , yo - A 0) denote Dirac's distribution, location of the load, its amplitude and 
excitation frequency, respectively. Once determined the transfer function (labeled as TF) 


corresponding to a given excitation frequency, its counterpart in the time domain can be obtained 
as the inverse Laplace transform £ : 

TF{x,y,t)= £- l {TF{x,y,s)}=^-F' + ~TF{x,y,s)e s 'ds. (12) 

In addition to the direct role in the determining of the response, the transfer function TF has then 
the role in determining the response to arbitrary time-dependent external excitations. 

The general procedure to identify the aeroelastic kernels of various order (bn), is to consider a 
general input in the form of Eq. (10) and to equate, for the generic term of n order, the 
coefficients of X,X 2 ••• X , l e {s ' +S2+ "' +s " ) ‘ . As an example, the first aeroelastic Volterra’s kernels 
that describes the linear system, obtained by neglecting the nonlinear terms in the aeroelastic 
governing equations, is obtained by considering the input load as L h {t')— X { e 1 (which in 
dimensionless form is expressed as l b {t)=(b/mUl)x i e' [ 'y, the response of the system is 
postulated in the form h{t)= H l {s l )X l e^ + h.o.t . Substituting h{t) and its derivatives in the 
governing equation of motion, one determines the coefficient of X,e . 

In a linear aeroelastic formulation, the system is completely characterized by a transfer function 
H { {sj) that contains the aerodynamic term as follow: 

#,(*,)= ( k hi +ms \ 2 +C >,1 S \ + s iP C La b U~C[- isfi/u^ ]+ j pc^y} . (13) 

Herein the Theodorsen's function C, connected with the Wagner's indicial function <p(t) via the 

Laplace's transform as C(- is)- (p{t)e' sr dr , has been included in the formulation. The terms 

underscored by the solid line correspond to the unsteady aerodynamic loads component 
(circulatory term), while the dotted line corresponds to the effect of the added mass. When the 
aerodynamic loads are neglected and for s = ico, this result coincides with that of the linear FRF, 
derived via the conventional Modal Analysis. Notice that, the condition s = ico correspond to 
zero initial conditions; the effect of the initial condition on the nonlinear aeroelastic response can 
be included remembering the complete form of the complex variable s, in which also the real 
part of the variable is included, namely, s = <7 + ico . 

For purely mechanical systems, in the frequency domain, the response via Volterra series has 
been carried out by several authors. In the present study an alternative procedure, based on the 
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multivariable kernel transforms techniques, referred to as Higher-order Transfer Functions 
(HTFs) is pursued. The two above mentioned approaches can be correlated each other, and this is 
shown also in this work. Assuming zero initial condition, the frequency response functions 
(FRFs) are obtained from the transfer functions TFs, by replacing the Laplace transform variable 

s with j(0 (Worden et al. 36 ). 

In the present nonlinear aeroelastic system, toward the estimation of higher order frequency 
response functions (HO FRFs) that are defined as the multi-dimensional Fourier Transform 
(MDFT) of the Volterra’s kernels, a sequence of transfer functions are employed. 

By the use of the linear frequency-response-function H,(s ,) , the behavior of the linear system is 

easily determined. It will be necessary to find a complete set of Vol terra kernel transforms 
H n (s,,s 2 ,---s„) for nonlinear systems and for this, in practice, we will use a convergent truncated 


series. 

Probing the system with a single harmonic yields only the information about the value of the 
transfer functions terms on the diagonal line of the plane j,,s 2 , in the Laplace transformed 
space, where 5 , = s 2 . However, in order to obtain information elsewhere in this space, one 
should use multi-frequency excitations. 

In the same way, the second order Volterra Kernel can be determined applying a load depending 
on two different frequencies expressed as: L b {t)= X x e« + X 2 ^' . In this case we can express the 


2 s*>t 


plunging response in the form: 

*(/)=/*,(*, )X i e 5 ' 1 +H { (s 2 )X 2 e^ +H 1 {s i ,s l )Xfe 2s '' +H 2 {s 2 ,s 2 )X 2 2 e : 

+ f/ 2 (5 l ,5 2 )X,X 2 e (l « +I ’- ) '+// 2 (^^ 1 )X 2 X ie ^ +s ' > +^. . (14) 

Substituting Eq. (14) in Eq. (8) and equating the terms containing X 1 X 2 e^ 1+S2 ^ the second 
order aeroelastic Volterra Kernel in the Laplace transformed space is obtained: 

H 2 {s [ ,s 2 )= -{w»i +k h2 )H i (j, )H , (j 2 )H X (5, + s 2 ) , (15) 


where: 

H l {s l +s l )=ifc M +{s l +s i fm + c H {s } + s 2 ) 

+ ( Sl +s 2 )pC La bUS [~ ifo + *2 >/^~ ]+ 7 pC La (*, + ) 2 b 2 ) ' , (16) 

is the first order Volterra Kernel in the Laplace transformed space at the frequency ro, +co 2 (that 
is obtained from Eq. (13) in which 5 , is replaced by 5, + s 2 ). Notice that the terms H 2 (s l ,s l ) and 
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H 2 {s ? , Sl ) can be determined from Eq. (15) replacing * 2 with 5, and vice versa, respectively. 
Following the same steps, applying the load L b (t)=X l e s ' 1 + X 2 e Sl ‘ + X } e 5,t , equating the terms 
in the form X,X 2 X/ s '" h ^ and remembering that: 

H l {s ] +s 2 +s 2 )=(k l +(s l +s 2 + sjm + c l {s l +s 2 +s 2 ) 

+ (s, + s 2 + s 3 )pC La bUS[- it. + *2 + *3 )b/Uj+\ pb 2 C La { S] +s 2 +s 2 f}\ (17) 

the expressions for the third order Volterra Kernel in the Laplace transformed space results: 

H 2 ($, , S 2 , S 2 ) = {H ] (s 3 1 (•*! {s 2 X^/.3 + C h3 S l S 2 S i ) 

+ 2H ,{s x ,s 2 \k h2 +c h2 {s t +s 2 )s 2 ))+2(H t (s 2 )H 2 (s J ,s i \k h2 + c h2 s 2 (s l + s 3 )) 

- H x (s 2 )H 2 (s 2 , ,, Xi k„ + c h2 s l (s 2 + * 3 ))))/ (//, fo + s 1+ s 3 )) .(18) 

Notice that the constants k h2 and c A2 multiply the whole expression of H 2 , and this term 

vanishes if the quadratic nonlinear term is absent in the aeroelastic governing equations. Herein 
one of the general properties of Volterra’ s series is recalled, namely that if all nonlinear terms in 
the equation of motion for the system consist of odd powers of x and y, then the associated 
Volterra series have no even-order kernels, and as a result it will possess no even-order TFs. It is 
also a general property of systems that all higher-order TFs can be expressed in terms of //, . 
The expressions are function of the system and can be obtained by using the harmonic probing 
algorithm. It clearly appears that the higher order of FRFs, defined from the Volterra series, are 
independent of the input to the system. 

C. Plunging-Pitching Airfoil Motion in an Incompressible Flow Field 

The governing aeroelastic system of an airfoil featuring plunging - twisting coupled motion, 
exposed to a harmonic time dependent external excitation is: 

mti + S a a + ^(c hi h' +k hj h‘)-L a =L b , ( 19 ) 

1=1 

S a h + I a d + )~ M “ = 0 • (20) 

i=l 

Considering, as usual, the blast load L b (r) as uniformly distributed in the chordwise direction, no 
moment contribution M b {t) is introduced in Eq. (20). 
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Following the steps adopted for 1-DOF, applying a load depending on one frequency L b -Xy , 
and expressing the plunging and pitching displacements in terms of transfer functions as. 

/ l (r)=X I // 1 ' , (5 l y , ' + X 1 2 // 2 A (^ 1 ,5 l > 2 ' , ' + X^H^{s„s l ,sy Sl ‘ , (21) 

+X^H“(s l , S y s '' +X^H^{s x ,s l ,s i y s \ ( 22 ) 

the relative kernels and the aeroelastic responses can be determined. 

The aeroelastic governing system including the blast pressure signatures can be expressed in the 
Laplace transformed space as: 

+ 5>( s ) + l.j 2 (|-aa)+-^sa = / ( ,(s)» (23) 
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Herein (•)= .£(■), and consequently <f = ^( 0 ) an d a = £(a(f))- 

Following the same steps, applying the loads L b {t)= X x e s ' 1 +X 2 e Sl ‘ and 

L b (t)=X x e s '’ + X 2 e s ' J + X 3 e h ', equating the terms in the forms X 1 X 2 e (j,+ ^ > and 

X l X 2 X 3 e < ' s ' +S2+s ^ t , the expressions for the second and third order Volterra Kernel in the Laplace 


transformed space can be obtained. 

D. Generalization to Multi Degrees of Freedoms Systems (M-DOFs) 

The method shown for 1-DOF and 2-DOF lifting surfaces can be extended to systems featuring 
multi degrees of freedoms in general, and to a 3 -D aircraft wing in particular. The method of 
deriving the n-f/i order nonlinear aeroelastic transfer functions is based upon the fact that when 
the aeroelastic system described by the response y(r) (expressed via Volterra series), is excited 


by a set of k unit amplitude exponentials at the arbitrary frequencies s h s 2 , ■■ 
contain exponential components of the form: 
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The presence of nonlinearities causes harmonic excitations and sums of harmonics to appear in 
the response of the aeroelastic system. Due to the nonlinear formulation, different frequencies 
can be expected as well. 

From the energetic point of view, we can observe that //, (s, ) produces a single frequency output , 
in response to the simple input c s ' . However, because the system is nonlinear, takes 

into account the terms that produce an output energy corresponding to the sum of frequencies 
60, +W-, , or in other words to the input e (l ' +l2V . Similarly, the third order nonlinear aeroelastic 
kernel, will inject a mix of three input frequencies into the total system output (see Figs. 6-8). 
This is the great advantage of this methodology over the other approaches based on the tirst 
order frequency response function. In contrast to these methodologies Volterra’s series approach 
is able to capture the transfer of energy between frequencies, that is typical for nonlinear 
systems. 

V. Results and Discussion 

To assess the versatility and provide a validation of this methodology, a comparison of the 
predictions of the aeroelastic response of nonlinear 2-D lifting surface using three 
approximations are shown in Figs. 2 and 3. The excellent agreement of the predictions, assess 
both the accuracy of the aeroelastic model and also the power of the methodology based on the 
Volterra series and indicial function approach. The first, the second and the third approximations 
of the aeroelastic response to the two loads 1 -COSINE gust load and triangular blast load are 
plotted for different parameters, together with the “exact” response of the aeroelastic system as 
obtained through digital-computer solution of the nonlinear aeroelastic governing equations. 
Both figures reveal the rapid convergence of the approximation. The same approach has been 
applied to a nonlinear time-varying system represented in Ref. 37 in which a transient response 
analysis of a continuous system has been addressed via functional techniques and 
multidimensional Laplace transformation. The impulse response of the system, represented by 
the differential equation: 

+ Ac(t)+ km(t)c{t)+ £ c 2 (t)= R S{t), (26) 

dt 

evaluated with the present analysis coincides with that shown in the Ref. 37 in which the 
parameters in use are: A = 1,/? = i,k = 1, 10, £ = l,m(/) = sin(60 o t),60 o = 2jt,20k . Figs. 4 and 5 


16 


show the excellent agreement of these two approaches. In Figs. 5 the phase-space of the two 
responses have been supplied. 

The results provided in these figures constitute a strong test of the speed of convergence of the 
present method 3 ^ ,3 ^. The coefficients of the nonlinear elements are not small, nor is the period of -• 
the time-varying parameter long compared with the natural time constant of the system. For the 
2-D lifting surface encompassing pure plunging, the first three aeroelastic kernels in magnitude 
and phase are depicted in Fig. 6 as a function of the frequency, considering the representation 
along the diagonal of the plane (o v co 2 , i.e. for co = <u, = co 2 = <o 3 • As is clearly seen, a reduced 
influence on the response of the third kernel is experienced. 

In Figs. 7 the Volterra’s kernels for the lifting surface featuring plunging - pitching coupled 
motions are depicted. Also in this case the plots include the magnitude and phase for the kernels 
in plunging H- and pitching , in which i identifies the order of the kernel. In Figure 8, 3-D 
plots of the magnitude and phase of the second order kernel vs. the two frequencies <y, and <u 2 
are provided. The contour plots reveal the symmetry of this kernel respect to the diagonal 
represented by w, = a ) 2 . The third order Volterra kernel for the case in which 0) y = cu, are 

depicted in a 3-D plot in Figs. 9. 

A. Response in Time and Frequency Domains 

Determination of the subcritical aeroelastic response to any time-dependent externally applied 
load is useful in the design of wing structures and of the associated feedback control systems. In 
some circumstances of the nonlinear analysis we are only interested in the special case implying 
T| - x 2 - ■ ■ • = T„ = T I6 ' 39 . This case can be represented as: 

£ ( T )= ^ ( T i > r 2 ' x n )| ri=r|= ... =Tii=T • (27) 

This function g(i) has a corresponding Laplace transform G(s) (also called associated 
transform) in the single-dimensional Laplace transform space: G{s) = £ fe(r)]. The response in 
time can be obtained from H{s v s 2 ,---,s n ) by determining G{s) first and evaluating the single 
dimensional inverse Laplace transform g{ t). This approach is called association of variable 
The nonlinear aeroelastic response in the time domain is depicted in Figs. 10 for a 2-D lifting 
surface featuring the plunging degree of freedom. In this figure the first plot represents the linear 
impulse response that corresponds to the convolution integral for the linear analysis. The other 
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three plots represent the components of the response due to the second and the third order kernels 
and the total response as a combination of the three partial responses. The influence of the 
stiffness and of the damping coefficients on the response are displayed in Figs. 11 and 12. An 
increase of the nonlinear damping or of the stiffness coefficients contributes to the decrease of 
the magnitude of the kernels and consequently, of the response amplitude. This shows that the 
nonlinearities in the stiffness and damping play a beneficial role on the subcritical aeroelastic 
response 18 . Figure 13 highlights the effect of the speed parameter V (sU/bco a ) on the lifting 

surfaces subjected to sonic-boom pressure signature as shown in the inset. Herein T p denotes the 
positive phase duration of the pulse measured from the time of impact of the structure, v denotes 
the shock pulse length factor. For r — 2 the N-shaped pulse degenerates into a symmetric sonic- 
boom pulse, in the sense that its positive phase has the same characteristics as its negative one, 
and for r = 1 a triangular pulse that corresponds to an explosive pulse is obtained. It becomes 
apparent that the amplitude of the response time-history (that have been evaluated for practical 
use with three kernels) increases with the increase of V . Moreover, in a certain range of speeds, 
as time unfolds, a decay of the amplitude is experienced, which reflects the fact that in this case 
the subcritical response is involved. However, for the dimensionless speed parameter V greater 
then the flutter speed (this one was determined using the linearized aeroelastic system), the 
response becomes unbounded implying that the occurrence of the flutter instability is impending. 
Also in this case the nonlinear stiffness and damping coefficients play a beneficial role on the 
subcritical aeroelastic response. 

VI. Conclusions 

In this work several issues related to the approach of the nonlinear aeroelastic response via 
Volterra’s series approach have been presented. Following the same approach presented here, the 
character of the instability boundary, i.e. benign or catastrophic can also be addressed. It was 
also shown that the method based on Volterra series opens large opportunities toward 
approaching in an unified and efficient way problems of nonlinear aeroelastic response and 
flutter instability. Moreover, this approach can be extended as to include also active control 
capabilities. Comparisons of results carried out via Volterra series in conjunction with indicial 
functions approach and classical approach have been provided in Fig. 1 for the linearized model. 
For the full nonlinear model, such validations have been displayed in Figs. 2-5, in which a 
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comparison with the exact digital evaluation of the nonlinear aeroelastic response has been 
supplied. In addition, it is worth noting that the aerodynamic indicial functions (for 
incompressible/compressible flow fields) considered in conjunction with Volterra's series 
approach can be used as a powerful analytical tool for developing unsteady aerodynamic models , 
and a unified nonlinear aeroelastic model. To the best of the authors’ knowledge, with the 
exception of this paper, the problem of the aeroelastic response of lifting surfaces to external 
pulses via Volterra’s series and indicial function approach was not yet addressed in the literature. 
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Figure Captions 

Fig. 1 Aeroelastic response to Dirac delta impulse, as represented in inset. Comparison of 
response prediction based on the first Volterra kernel and the exact solution 
(m = l;c M =0.1;**, =10 2 £ = l;p= 0.125; U_ =0.4U F ,C La = 2 n). 

Fig. 2 Convergence study involving the first three kernels and comparison with the “exact” 
nonlinear aeroelastic response to a 1 -COSINE gust pulse, as shown in the inset; parameters as in 
Fig. 1. 

Fig. 3 Convergence study involving the first three kernels and comparison with the “exact” 
nonlinear aeroelastic response to a triangular blast load, as shown in the inset; parameters as in 
Fig. 1. 

Fig. 4 Convergence study involving the first three kernels and comparison with the “exact” 
nonlinear aeroelastic response of the impulse response described in Eq. (26) (a = \-,R = 1;* = 1; 

£ = l;£o 0 = 2 k )- 

Fig. 5 Convergence study involving the first three kernels and comparison with the “exact” 
nonlinear aeroelastic response of the impulse response described in Eq. (26) (phase-space 
representation) (A = 1;/? = l;Jk = 10;£ = 1; fO 0 = 20?r). 

Fig. 6 Comparison of the first three aeroelastic kernels for pure plunging motion 
(m = l;c„ = 10;*„, =10 4 ;c w =10;* M =10 7 ;c„ = 10;**, =10 ,0 ;6 = l;p =0.125; U„ =0.4t/ f \C La = In). 

Representation for s i = s 2 = s 3 , i.e. £0, = 0) 2 =co i . 

Fig. 7 First two aeroelastic kernels for plunging - pitching coupled motion (m = l;a = -0.2; c hl = 10; 
(ifc A1 = 10 4 ;c al -10;^, =10 4 ;c w =10;* A2 =10 7 ;c a2 =10;* o2 =10 7 ;fc = l;p = 0.125; [/„ =0.4t/ f -,C La = 2n). 
Fig. 8 3-D and contour plots of second order aeroelastic kernel; parameters as in Fig. 6. 

Fig. 9 3-D and contour plots of third order aeroelastic kernel; parameters as in Fig. 6. 

Fig. 10 Time-history of the nonlinear aeroelastic response (U_ =0.4 V F \C La =2n\m = V,c hi =10; 

k M =10 4 ;c A2 =10; k h2 =10 7 ;c„ =10;* M =10 ,0 ;fc = l;p= 0.125). 

Fig. 11 Influence of the nonlinear stiffness coefficient k hi on the nonlinear aeroelastic response; 
parameters as in Fig. 7. 

Fig. 12 Influence of the nonlinear damping coefficient c*, on the nonlinear aeroelastic response; 
parameters as in Fig. 7. 

Fig. 13 Influence of the flight speed on the nonlinear aeroelastic response to a sonic-boom 
(r p = 15sec;r = 2 ), as shown in the inset, evaluated with three kernels; parameters as in Fig. 7. 
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